A Study on a Radio Source Location Estimation System Using High Altitude Platform Stations (HAPS)

Currently, there is a system in Japan to detect illegal radio transmitting sources, known as the DEURAS system. Even though crackdowns on illegal radio stations are conducted on a regular basis every year, the number of illegal emission cases still tends to increase, as ordinary citizens are now able to handle advanced wireless communication technologies, e.g., via software-defined radio. However, the current surveillance system may not be able to accurately detect the source in areas where large buildings are densely packed, such as urban areas, due to the effects of reflected waves. Therefore, in this study, we proposed a system for estimating the location of the source of transmission using a high-flying unmanned aerial vehicle called HAPS. The simulation results using numerical analysis software show that the proposed system can estimate the location of the source over a wider area and with higher accuracy than conventional monitoring systems.


Introduction
In recent years, mobile communication technology has made great progress, and is used not only in public places but also in our everyday life.For example, mobile communication technology is indispensable in the field of IoT, which is used in smart phones and electrical appliances that everyone has today, and in the field of automatic driving.In fact, the penetration rate of mobile devices is increasing every year [1].The number of wireless stations has been increasing along with this growth, and as of the end of September 2022, a total of nearly 3 million stations were in operation nationwide in Japan, as shown in Figure 1 [2].As time has progressed, technologies such as cell phones and IoT have been developed, and mobile communication devices have become more familiar to us, the general public.However, perhaps due to these developments, the number of illegal radio stations (which use power in excess of the limit or illegal frequencies without a license) has been increasing in recent years, as shown in Figure 2 [3].As time has progressed, technologies such as cell phones and IoT have been developed, and mobile communication devices have become more familiar to us, the general public.However, perhaps due to these developments, the number of illegal radio stations (which use power in excess of the limit or illegal frequencies without a license) has been increasing in recent years, as shown in Figure 2 [3].As time has progressed, technologies such as cell phones and IoT have been developed, and mobile communication devices have become more familiar to us, the general public.However, perhaps due to these developments, the number of illegal radio stations (which use power in excess of the limit or illegal frequencies without a license) has been increasing in recent years, as shown in Figure 2 [3].Today, these illegal radio stations are searched in various ways, and one of them is called a radio wave monitoring system.This system monitors radio waves received at sensor stations and estimates the direction of the radio wave source.Based on this, the location of illegal radio stations is identified by remotely operating sensor stations installed at various locations and sensor stations mounted on vehicles, from a center station installed at each general communications station [3].However, in urban areas where many buildings are densely built or where there is a lot of vegetation, these systems may not work well because there may be no line of sight from the sensor station to the illegal radio station, or many reflected waves may be generated.In recent years, research has been conducted using unmanned aerial vehicles (UAVs) such as drones to gain some altitude in order to improve the location estimation accuracy of the source of radio waves by making Today, these illegal radio stations are searched in various ways, and one of them is called a radio wave monitoring system.This system monitors radio waves received at sensor stations and estimates the direction of the radio wave source.Based on this, the location of illegal radio stations is identified by remotely operating sensor stations installed at various locations and sensor stations mounted on vehicles, from a center station installed at each general communications station [3].However, in urban areas where many buildings are densely built or where there is a lot of vegetation, these systems may not work well because there may be no line of sight from the sensor station to the illegal radio station, or many reflected waves may be generated.In recent years, research has been conducted using unmanned aerial vehicles (UAVs) such as drones to gain some altitude in order to improve the location estimation accuracy of the source of radio waves by making it easier to see the source of radio transmission [4].However, the service area of this system is limited by the height of the UAVs and actual implementation has not yet been achieved.As for radio monitoring technology, advanced spectrum sensing techniques are still being developed to enhance the detection of illegal radio transmissions in harsh environments [5].Many other studies on radio monitoring technology [6,7] still exist today.
To address this issue, this paper focuses on the use of the High Altitude Platform System (HAPS) which is a generic term for a system that provides communication services over a wide area by operating an unmanned aircraft, such as an aircraft flown in the stratosphere, like a communication base station [8,9].As explained here, the HAPS mainly operates unmanned aircraft as base stations, and papers such as "Research and Development of HAPS Mobile Communication System for Rapid Disaster Recovery" [10] have been published.The HAPS has the advantage of being able to deploy unmanned base stations regardless of ground traffic conditions or environments such as the sea.So far, research has been conducted with the idea of operating unmanned aircraft as the main base station.However, if the HAPS is operated as an illegal radio surveillance system, it will be possible to monitor a wider area with less influence from reflected waves and at a lower cost than conventional surveillance systems.If the HAPS can also be operated as a base station in times of disaster, it will be effective in estimating the location of mobile terminals in times of disaster.This paper examines the operation of the HAPS as an illegal radio monitoring system and a location estimation system for mobile terminals during disasters, and demonstrates its effectiveness mainly through numerical simulations using MATLAB [11].This paper aims to evaluate the performance of a location estimation system using the HAPS for illegal wireless stations and mobile terminals during disasters that occur in various environments and distances through numerical simulations.

Methods for Estimating the Location of Radio Transmission Sources
There are two types of systems for estimating the location of radio wave sources such as illegal radio waves and mobile terminals: geometric methods that use various numerical data observed from actual radio wave sources to estimate the location, and statistical methods that estimate the location by comparing the actual observation data with a database obtained using a virtual transmitter in advance.This section describes the geometrical methods currently used to estimate the location of a radio transmitter.In Japan, the DEURAS system [3] is currently used to monitor illegal radio stations.The DEURAS system uses phase differences of the received signals at different antennas to estimate the direction of incoming radio waves.This section briefly explains the methods based on estimation of angle of arrival.

Estimation Method Using Phase Difference of Received Signals
This section revisits a well-known method for estimating the location of a radio wave source.The direction of arrival of a radio wave is determined by using the phase difference between received signals at multiple elements based on the same principle as that of interferometry.For simplicity, we will explain the estimation of the direction of arrival when a plane wave is incident on two array antennas, as shown in Figure 3.When considering such a situation, the phase difference ∆ϕ between the received signals of arrays 1 and 2 can be expressed, as in Equation ( 1) [12].If the complex received signal series of each antenna are  and  , respectively, the phase difference ∆ can be expressed as in Equation (2) [12].

∆𝜙 = tan
Im    * Re    * where E[•] denotes an ensemble average operation,  is the wavelength of the carrier wave,  is the antenna spacing, Re • donotes the real part, and Im • denotes imaginary part.When the antenna spacing is 1/2 of the carrier wave wavelength,  = sin ∆ from Equation (1), and the direction of arrival  can be obtained in the range −  by determining the phase difference ∆.In a previous study [13], in order to improve the estimation accuracy, the direction of arrival was obtained by using several pairs of two-element arrays to obtain (in advance) the complex reception response  () If the complex received signal series of each antenna are x 1 and x 2 , respectively, the phase difference ∆ϕ can be expressed as in Equation (2) [12].
where part.When the antenna spacing is 1/2 of the carrier wave wavelength, θ = sin −1 ∆ϕ π from Equation (1), and the direction of arrival θ can be obtained in the range − π 2 < θ < π 2 by determining the phase difference ∆ϕ.In a previous study [13], in order to improve the estimation accuracy, the direction of arrival was obtained by using several pairs of two-element arrays to obtain (in advance) the complex reception response a i (θ) to the arriving wave from the angle θ of element i and introducing the evaluation function shown in Equation (3) [12].When P(θ) reaches its maximum, θ is consistent with the direction of arrival.
The estimated direction of arrival is then used to estimate the location of the radio source by triangulation, as shown in Figure 4.
If the complex received signal series of each antenna are  and  , respectively, the phase difference ∆ can be expressed as in Equation ( 2) [12].

∆𝜙 = tan
Im    * Re    * where E[•] denotes an ensemble average operation,  is the wavelength of the carrier wave,  is the antenna spacing, Re • donotes the real part, and Im • denotes imaginary part.When the antenna spacing is 1/2 of the carrier wave wavelength,  = sin ∆ from Equation (1), and the direction of arrival  can be obtained in the range −  by determining the phase difference ∆.In a previous study [13], in order to improve the estimation accuracy, the direction of arrival was obtained by using several pairs of two-element arrays to obtain (in advance) the complex reception response  () to the arriving wave from the angle  of element i and introducing the evaluation function shown in Equation (3) [12].When () reaches its maximum,  is consistent with the direction of arrival.
The estimated direction of arrival is then used to estimate the location of the radio source by triangulation, as shown in Figure 4.This method provides highly accurate estimation when there is a line of sight (LoS) from the antenna to the source of the radio wave.However, when estimating the direction of arrival in urban areas, multiple waves are received due to reflected waves.The problem with this method is that the direction of arrival cannot be estimated properly at that time.

Estimation Method Using the Beamformer Method
The Beamformer method is a basic method for estimating the direction of arrival.It estimates the direction of a radio source by scanning the main lobe of a directional antenna with high gain in all directions, as shown in Figure 5, outputting the output voltage change of the array at each angle, and obtaining the maximum output value.This method provides highly accurate estimation when there is a line of sight (LoS) from the antenna to the source of the radio wave.However, when estimating the direction of arrival in urban areas, multiple waves are received due to reflected waves.The problem with this method is that the direction of arrival cannot be estimated properly at that time.

Estimation Method Using the Beamformer Method
The Beamformer method is a basic method for estimating the direction of arrival.It estimates the direction of a radio source by scanning the main lobe of a directional antenna with high gain in all directions, as shown in Figure 5, outputting the output voltage change of the array at each angle, and obtaining the maximum output value.For simplicity, the estimation of the direction of arrival using the Beamformer method is described in detail using a K-element linear array antenna, as shown in Figure 6 [12].For simplicity, the estimation of the direction of arrival using the Beamformer method is described in detail using a K-element linear array antenna, as shown in Figure 6 [12].For simplicity, the estimation of the direction of arrival using the Beamformer method is described in detail using a K-element linear array antenna, as shown in Figure 6 [12].In this case, the output of the array is () can be expressed as follows [12].In this case, the output of the array is y(t) can be expressed as follows [12].
(•) H denotes the transpose operation, where the rows and columns of a matrix are swapped.
(•) T denotes the Hermitian transpose operation, where the matrix is first transposed and then each element is replaced by its complex conjugate.The output voltage P out is provided as follows [12].
Assuming that L waves are incident on the incoming wave and that the respective signal waveforms and angles of arrival are F l (t) and θ l (l = 1, 2, . . ., L), the respective direction vectors V l can be expressed as follows [12].
Also, d k (k = 1, 2, . . ., K) is the distance from the reference array to the position of the k-th element.At this time, the input vector is expressed as follows [12].
The matrix S in Equation ( 16), which represents the correlation between arriving waves, is called the signal correlation matrix and is a diagonal matrix with P l (l = 1, 2, . . ., L) (each input power) in the diagonal components, if all arriving waves are uncorrelated with each other.In order to direct the main lobe of the array antenna to angle θ using the above equation and others, the weight w k expressed in Equation ( 17) [12] can be multiplied to each input signal based on the common-phase condition (the condition to align the phases to be in-phase).
The angle θ is varied from −π/2 to π/2 to find the peak of the output voltage of the array.The output voltage of the array is expressed by Equation ( 18), and the angular distribution by the Beamformer method is obtained by normalizing Equation ( 18) [12].
The Beamformer method estimates the direction of arrival and the input power of the arriving wave from the position of the peak of the evaluation function, as shown in Equation ( 18) using the input correlation matrix R xx and the weight expressed in Equation (17).Using the estimated direction of arrival, the radio wave source is estimated by triangulation, as shown in Figure 4.However, in the case of multiple incoming radio waves, if the main lobe of the antenna pattern shown in Figure 5 faces a certain radio wave while the side lobes (in the direction of relatively high gain) coincide with the direction of other radio waves, the output voltage will include not only one wave but also other waves, and it is not guaranteed that the evaluation function may always have a peak in the desired direction of arrival.In addition, when radio waves with almost similar arrival angles are incident, it may not be possible to distinguish the two waves when scanning with the main lobe having a wide angle of arrival, as shown in Figure 5.In such cases, the angular resolution of location estimation by the Beamformer method becomes low, and the peak value may not accurately represent the power of the signal.
In fact, as shown in [14] which studied the direction of arrival estimation using the Beamformer method in a multipath environment, the estimation accuracy is degraded depending on the number of arriving waves.There is also the Capon method, which directs the main lobe in one direction while at the same time attempting to reduce the output contribution from other directions.However, all of these methods use the main lobe of the antenna pattern, and the angular resolution of the estimation is limited.Therefore, it is difficult to apply the Beamformer method to the scenarios of ground emitters in urban areas where many radio waves arrive.

Estimation Method Using the MUSIC Method
This section describes in detail the MUSIC method, which uses eigenvalues and eigenvectors to estimate the direction of arrival at a higher angle resolution than the Beamformer method [12].For the sake of simplicity, let assume a K-element linear array antenna is used, as shown in Figure 6.In this case, the input vector X(t) of each antenna is the same as in Equation (12).
The correlation matrix of the input vectors is then as in Equation ( 15).
where the signal correlation matrix S can be expressed as follows.
In the absence of thermal noise, if the arriving waves are uncorrelated, S becomes a diagonal matrix and its rank is L. If the direction matrix A also differs from the direction of arrival, its column vectors become independent, and its rank becomes L. Therefore, in this case, the correlation matrix R xx is a non-negative definite Hermitian matrix of rank L. Denoting the eigenvalues of this matrix by µ i (i = 1, 2, . . ., K) and the corresponding eigenvectors by e i (1, 2, . . . ,K) [12], and each eigenvalue can be expressed as a real number [12].
The eigenvectors are orthogonal to each other, as follows [12].
In the presence of thermal noise, the thermal noise power is added to the eigenvalues and the eigenvectors remain unchanged [12].
Then we have: Denoting the eigenvalues of the correlation matrix by: the number of arriving waves L can be obtained by finding the number of eigenvalues of the correlation matrix that have values larger than the noise power.Also, if we focus on the eigenvectors corresponding to the eigenvalues equal to the thermal noise power: based on Equation (20), we can derive the following property [12].
Also, from the fact that both A and S are full ranks.
In other words: From this equation, it can be seen that the eigenvectors corresponding to eigenvalues equal to the thermal noise power are orthogonal to the direction vectors of all arriving waves.Using this property, the evaluation function of the MUSIC method can be expressed as follows [12] with normalization.
The denominator part of Equation (32) shows a minimum value at θ, which is the direction vector (direction vector of the arriving wave) orthogonal to the eigenvector corresponding to the noise, as mentioned earlier.The MUSIC method, on the other hand, uses the noise portion (the portion where the power is small) to achieve high resolution estimation.However, when multiple radio waves are incident from a single source due to reflected waves, etc., there is a correlation between the incoming waves.In this case, the rank of the signal correlation matrix expressed in Equation ( 21) is not equal to the number of arriving waves, and the direction of arrival cannot be accurately estimated.
It is known that all of the three location estimation methods using direction of arrival estimation introduced in this section have the problem that the estimation accuracy is degraded by many reflected waves when a ground antenna is used.

Numerical Analysis and Results
Currently, location estimation using the HAPS is not yet available and is considered to be operated as a base station.Because of this situation, no demonstration experiment of location estimation using the HAPS has been realized.In this section, two methods of direction of arrival estimation will be applied for the HAPS for the location estimation of ground emitters.Since the phase difference method is just a primitive version of the Beamformer one, only the Beamformer and MUSIC methods are introduced in this section for evaluation of the localization accuracy.

Comparison of LoS Rate between the HAPS and UAVs
In Section 2, we discussed various location estimation methods.However, the accuracy of such methods varies greatly depending on whether the environment between the receiving antenna and the radio source is a LoS environment or an NLoS environment.In Sensors 2024, 24, 5803 9 of 22 this section, we perform a ray tracing simulation using the SBR (Shooting and Bouncing Rays method) in MATLAB to compare the LoS rate between the case using the HAPS and the case using a UAV.In this analysis, the case in which a direct wave (a radio wave that is not reflected by any obstacle) reaches the receiving antenna was considered as the LoS environment and the case in which direct waves do not reach the receiving antenna considered an NLoS environment.

Simulation Model
Simulations were performed using a pre-existing MATLAB 3D model environment that mimics the town of Chicago, shown in Figure 7.An overall view of the simulation is shown in Figure 8.The simulation was conducted assuming that the HAPS was 2000 m above the ground and the UAV was 200 m above the ground.For simplicity, we will use the simulation of the LoS rate over the town to illustrate the simulation.First, random radio sources are generated in the town, as shown in Figure 9. Here, the height of each antenna was set to 10 m.An overall view of the simulation is shown in Figure 8.The simulation was conducted assuming that the HAPS was 2000 m above the ground and the UAV was 200 m above the ground.For simplicity, we will use the simulation of the LoS rate over the town to illustrate the simulation.First, random radio sources are generated in the town, as shown in Figure 9. Here, the height of each antenna was set to 10 m.For simplicity, we will use the simulation of the LoS rate over the town to illustrate the simulation.First, random radio sources are generated in the town, as shown in Figure 9.
Here, the height of each antenna was set to 10 m.Next, as shown in Figure 10, a receiving antenna (in this case, a UAV) is generated over the town, the number of direct waves is calculated using a ray tracing simulation, and the LoS rate is calculated by calculating the ratio of the number of direct waves to the number of radio wave sources.In this simulation, the number of radio wave sources is set to 1000.The same simulation is performed for the HAPS.Next, as shown in Figure 10, a receiving antenna (in this case, a UAV) is generated over the town, the number of direct waves is calculated using a ray tracing simulation, and the LoS rate is calculated by calculating the ratio of the number of direct waves to the number of radio wave sources.In this simulation, the number of radio wave sources is set to 1000.The same simulation is performed for the HAPS.Next, as shown in Figure 10, a receiving antenna (in this case, a UAV) is generated over the town, the number of direct waves is calculated using a ray tracing simulation, and the LoS rate is calculated by calculating the ratio of the number of direct waves to the number of radio wave sources.In this simulation, the number of radio wave sources is set to 1000.The same simulation is performed for the HAPS.When the simulation is performed away from the center, receivers are generated in several directions, as shown in Figure 11, and the average value of the LoS ratio obtained from each receiver is used as the LoS ratio.When the simulation is performed away from the center, receivers are generated in several directions, as shown in Figure 11, and the average value of the LoS ratio obtained from each receiver is used as the LoS ratio.

Results of LoS Rate Measurement Simulation
The results of the simulation conducted above are shown in Figure 12.The LoS ratio of the HAPS was about 100% when the center of the urban area (where the radio source was generated) was directly below the HAPS, and about 60% at a distance of 20 km (elevation angle 45°).On the other hand, when a UAV was used, the LoS ratio was about 30% even when the center of the urban area was directly below the UAV, indicating that the majority of the radio sources were out of sight.Therefore, a specific altitude for the HAPS is considered necessary for position estimation by angle estimation.

Results of LoS Rate Measurement Simulation
The results of the simulation conducted above are shown in Figure 12.The LoS ratio of the HAPS was about 100% when the center of the urban area (where the radio source was generated) was directly below the HAPS, and about 60% at a distance of 20 km (elevation angle 45 • ).On the other hand, when a UAV was used, the LoS ratio was about 30% even when the center of the urban area was directly below the UAV, indicating that the majority of the radio sources were out of sight.Therefore, a specific altitude for the HAPS is considered necessary for position estimation by angle estimation.

Location Estimation Simulation
In order to realize a simulation of radio source location estimation using the HAPS, we used MATLAB to estimate the location of radio sources.Several location estimation methods were discussed in Section 2, but considering the wide search area and required

Location Estimation Simulation
In order to realize a simulation of radio source location estimation using the HAPS, we used MATLAB to estimate the location of radio sources.Several location estimation methods were discussed in Section 2, but considering the wide search area and required estimation accuracy, the Beamformer method and the MUSIC method were used in our simulation.As a comparison, a case in which location estimation was performed at the altitude of the UAV was also simulated.In our analysis, the location estimation simulation was performed for a non-mobile illegal radio source located outdoors in an urban area, as shown in Figure 7.
The simulation process is shown in Figure 13.(1) A radio source is randomly generated, (2) a ray tracing simulation is performed to obtain information about arriving radio (direction of arrival, number of reflections, propagation distance, etc.), and (3) the received signal is generated from the aforementioned information.The direction of arrival is estimated by using each estimation method applied on the received data generated in (3).Then, location estimation is performed using the triangularization method.By comparing the randomly generated coordinates of the radio source with the estimated coordinates, the location estimation of the radio source using the HAPS and other methods is evaluated.Figure 13 shows an overview of the simulation.
x FOR PEER REVIEW 14 of 24

Receiving and Transmitting Antenna Configuration
Table 1 shows the antenna specifications for estimating the location of a radio wave source using ray tracing simulation.Hereafter, the transmitting antenna is denoted as Tx and the receiving antenna as Rx.

Receiving and Transmitting Antenna Configuration
Table 1 shows the antenna specifications for estimating the location of a radio wave source using ray tracing simulation.Hereafter, the transmitting antenna is denoted as Tx and the receiving antenna as Rx.The antenna was modeled after the cylinder antenna used in the HAPS as a base station for the simulation.The actual prototype cylinder antenna can be found in [15].The specifications are shown in Table 2.In this paper, the cylinder antenna consists of N = 36 elements, as depicted in Figure 14, where each element was assumed to be omni-directional and the maximum gain was set to 0 dBi in the simulation.In this paper, the cylinder antenna consists of  = 36 elements, as depicted in Figure 14, where each element was assumed to be omni-directional and the maximum gain was set to 0dBi in the simulation.

Radio Propagation Model
This section describes the radio propagation model used in the simulation of radio source location estimation using the HAPS.In this simulation, statistical data were not used in the propagation model because a ray tracing simulation was conducted to determine specific radio waves.The free propagation loss was modeled using Friis formula, as follows [16].

Radio Propagation Model
This section describes the radio propagation model used in the simulation of radio source location estimation using the HAPS.In this simulation, statistical data were not used in the propagation model because a ray tracing simulation was conducted to determine Sensors 2024, 24, 5803 14 of 22 specific radio waves.The free propagation loss was modeled using Friis formula, as follows [16].
where d and λ are the distance between the transceiver pairs and the wavelength of the considering radio wave, respectively.The overall path loss can be expressed by Equation ( 36), where the additional attenuation L re f due to earth and building reflections was calculated using the Fresnel equation.The permittivity and conductivity of the surface materials used in these calculations were taken from the ITU-R study [17].Assuming that both the ground surface and the building are concrete, the relative permittivity and conductivity were calculated as 5.31 and 0.0548 [S/m], respectively.

Signal Processing System
This section describes the signal processing performed in the location estimation simulation.First, the analytical signal input to the cylinder antenna is obtained by generating a single radio source in a 3D model and considering the case where L radio waves are incident from the source.The propagation loss PL l [dB] (l = 1, 2, . . ., L) is obtained from Equation (37) using the propagation distance and other data obtained from the ray tracing simulation, and the input power is calculated by following formula.
Using this value, the propagation time delay against the direct wave τ l (l = 1, 2, . . ., L), (τ 1 = 0) and the times of reflection m l (l = 1, 2, . . ., L), the complex waveform of the l-th wave can be written by the following formula.
The array matrix for the cylinder antennas is then expressed by the following equation [12].
where r k (k = 1, 2, . . ., N) is the position vector of each element, and θ and ϕ are the elevation and azimuth angles to the reference point, respectively.The thermal noise power P N generated at the receiving antenna is given by following equation [18].
where k, T, B, NF are Boltzmann's constant, absolute temperature of the noise source, bandwidth, and noise figure, respectively.The absolute temperature in the HAPS environment was assumed to be 216.5 [K] and the bandwidth was 10 [MHz].The noise figure was calculated as 3 dB [19].As a result, the thermal noise power P N was about −102 [dBm] in the HAPS environment while in the UAV environment, the thermal noise power was slightly higher since the absolute temperature at the UAV altitude is about 288 [K].Using Equations ( 38), ( 39) and (42), the input vector X(t) to each antenna is expressed by the following equation.
In the above equation, N(t) is the thermal noise vector, whose components are independent complex Gaussian processes with mean zero and variance (power) of P N .The input vector, obtained by substituting appropriate values for t in Equation ( 43) for different time samples, was used as the input complex signal for the cylinder array.The number of samples was set to 100 in this paper.By applying the aforementioned Beamformer and MUSIC methods to the obtained input vector X(t), the direction of arrival from the radio source was estimated and used for the location estimation method described in the next section.The direction of arrival with the highest power obtained by the Beamformer and MUSIC methods was used as the estimated direction of arrival.

Location Estimation Method
This section describes a method for obtaining the estimated coordinates of a radio source from the direction of arrival obtained by the MUSIC and Beamformer methods.In this simulation, the direction of arrival of radio waves is estimated from three points as shown in Figure 15.The latitude and longitude of each HAPS observation point are converted to the Cartesian coordinate system, and a direction line is created in the x-y coordinate system from the direction of arrival.Figure 16 shows two localization methods from the estimated angular information.In the left-hand side case, all angular information is available from the three sensors, and the center of gravity of the intersections of each directional line is considered as the estimated coordinates of the emitter.In the right-hand side case, the angular information from one of three sensors is missing, so the estimated coordinates are the intersection points of the direction lines of the remaining sensors where angular information is available.For all other cases, the location of the emitter is considered to be unidentified (undetected).The latitude and longitude of each HAPS observation point are converted to the Cartesian coordinate system, and a direction line is created in the x-y coordinate system from the direction of arrival.Figure 16 shows two localization methods from the estimated angular information.In the left-hand side case, all angular information is available from the three sensors, and the center of gravity of the intersections of each directional line is considered as the estimated coordinates of the emitter.In the right-hand side case, the angular information from one of three sensors is missing, so the estimated coordinates are the intersection points of the direction lines of the remaining sensors where angular information is available.For all other cases, the location of the emitter is considered to be unidentified (undetected).gular information.In the left-hand side case, all angular information is available from the three sensors, and the center of gravity of the intersections of each directional line is considered as the estimated coordinates of the emitter.In the right-hand side case, the angular information from one of three sensors is missing, so the estimated coordinates are the intersection points of the direction lines of the remaining sensors where angular information is available.For all other cases, the location of the emitter is considered to be unidentified (undetected).

Results of Location Estimation Simulation
The 3D city model shown in Figure 7 was used to randomly generate 1000 radio wave sources, and the location of each source was estimated using the method described in the previous section.The horizontal distance between the center of the urban area with illegal radio sources and the hovering sensors was varied with different scenarios, as seen in our evaluation results below.Three scenarios of sensor altitudes, i.e., HAPS, UAV, and DEURAS, are also separately evaluated in this section. HAPS altitude

Results of Location Estimation Simulation
The 3D city model shown in Figure 7 was used to randomly generate 1000 radio wave sources, and the location of each source was estimated using the method described in the previous section.The horizontal distance between the center of the urban area with illegal radio sources and the hovering sensors was varied with different scenarios, as seen in our evaluation results below.Three scenarios of sensor altitudes, i.e., HAPS, UAV, and DEURAS, are also separately evaluated in this section.

■ HAPS altitude
Table 3 shows the results of localization using the HAPS at three points 5 km away from the urban area.The Beamformer and MUSIC methods were used to estimate the direction of arrival.The estimation probability denotes the percentage that the emitter can be identified (with localization errors) via our proposed system, as mentioned in Figure 16.The stochastics of the estimation errors were derived only from emitters, which were identified by our proposed method.From Table 3, it is observed that there was no significant difference in the results between the Beamformer and MUSIC methods, and 90% of the estimated locations were within 30 m of the true coordinates.Table 4 shows the results of the estimation using the HAPS at three locations 10 km away from the urban area.Since the sensors are more distanced from the sources compared to the 5 km case, we can observe a slight degradation in terms of localization accuracies as compared to these in Table 3.Furthermore, Table 5 shows the results of the estimation using the HAPS at three locations 20 km away from the urban area.Similar to the two aforementioned scenarios, neither the Beamformer nor the MUSIC method showed a significant difference in results.Also, the accuracy of the estimation was poorer at greater distances compared to the cases at 5 km and 10 km.Furthermore, only about 83.6% of the radio transmitting sources in urban areas could be obtained.The CDF of localization errors for these scenarios is summarized in Figure 17. UAV altitude Similarly, a simulation was conducted assuming a UAV, with the height of the receiving antenna changed to 200 m.In addition, because the effect of reflections is likely to increase at the altitude of the UAV, the simulation was performed considering radio waves that can be reflected up to four times.Tables 6-8 show the results of the estimation using UAVs at three locations (5 km, 10km and 20km away from the urban area, respectively).
In Table 6, only 20% of the total number of radio sources could be identified due to the low LoS ratio, as described in Section 3.1.The localization estimation error of the radio sources that could be identified was also larger than that of the HAPS case.21.1 Table 7 shows the results of the estimation using UAVs at three locations 10 km away from the urban area.As in the 5 km case, the number of radio sources that could be identified was small.However, the estimation error and the CDF90% value was slightly improved against the 5 km case.This is thought to be due to the fact that the 10 km environment has smaller multipath components compared to the 5km case, such that angular estimation has higher accuracy.

■ UAV altitude
Similarly, a simulation was conducted assuming a UAV, with the height of the receiving antenna changed to 200 m.In addition, because the effect of reflections is likely to increase at the altitude of the UAV, the simulation was performed considering radio waves that can be reflected up to four times.Tables 6-8 show the results of the estimation using UAVs at three locations (5 km, 10 km and 20 km away from the urban area, respectively).In Table 6, only 20% of the total number of radio sources could be identified due to the low LoS ratio, as described in Section 3.1.The localization estimation error of the radio sources that could be identified was also larger than that of the HAPS case.
Table 7 shows the results of the estimation using UAVs at three locations 10 km away from the urban area.As in the 5 km case, the number of radio sources that could be identified was small.However, the estimation error and the CDF90% value was slightly improved against the 5 km case.This is thought to be due to the fact that the 10 km environment has smaller multipath components compared to the 5 km case, such that angular estimation has higher accuracy.
Table 8 shows the results of the estimation using UAVs at three locations 20 km away from the urban area.Both the estimated probability and localization accuracy are degraded in this case.It is thought to be due to the fact that the sensors are located too far away, such that the arrival waves after multiple reflections are too weak compared to the noise level.
The CDF of localization errors for these different scenarios is summarized in Figure 18.
Sensors 2024, 24, x FOR PEER REVIEW 20 of 24 Table 8 shows the results of the estimation using UAVs at three locations 20 km away from the urban area.Both the estimated probability and localization accuracy are degraded in this case.It is thought to be due to the fact that the sensors are located too far away, such that the arrival waves after multiple reflections are too weak compared to the noise level.The CDF of localization errors for these different scenarios is summarized in Figure 18. DEURAS Finally, a model simulating the DEURAS system was created by installing the receiving antennas on the rooftop of buildings with good visibility in an urban area, as shown in Figure 19.Other system parameters were the same as those for UAVs.■ DEURAS Finally, a model simulating the DEURAS system was created by installing the receiving antennas on the rooftop of buildings with good visibility in an urban area, as shown in Figure 19.Other system parameters were the same as those for UAVs.Table 9 shows the evaluation results of position estimation simulations using the Beamformer and MUSIC methods where DEURAS sensors are located, as shown in Figure 19.Compared to location estimation using UAVs, the estimation accuracy is still poorer, except for the fact that DEURAS can provide a higher estimated probability since the sensors are close to the sources.79.9 For a fair comparison, we conducted new simulations that the HAPS and UAVs were deployed at three points at a distance of 500 m from the center of the urban area (just right above the area similar to the case of the DEURAS system).
Table 10 shows the evaluation results of the estimation using the HAPS at three locations 500 m away from the urban area.Similarly, Table 11 shows the results of the evaluation when estimating using UAVs at three locations 500 m away from the urban area.Table 9 shows the evaluation results of position estimation simulations using the Beamformer and MUSIC methods where DEURAS sensors are located, as shown in Figure 19.Compared to location estimation using UAVs, the estimation accuracy is still poorer, except for the fact that DEURAS can provide a higher estimated probability since the sensors are close to the sources.For a fair comparison, we conducted new simulations that the HAPS and UAVs were deployed at three points at a distance of 500 m from the center of the urban area (just right above the area similar to the case of the DEURAS system).
Table 10 shows the evaluation results of the estimation using the HAPS at three locations 500 m away from the urban area.Similarly, Table 11 shows the results of the evaluation when estimating using UAVs at three locations 500 m away from the urban area.Figure 20 shows a summary of the CDF of localization errors of all scenarios.These results show that the HAPS shows the best performance in terms of both estimated probability and localization accuracy.It is owing to the fact that the HAPS has a higher LoS ratio to guarantee a higher estimated probability, and also direct waves are more likely to arrive at the HAPS rather than multiple reflection paths that increase the angular estimation accuracy.Figure 20 shows a summary of the CDF of localization errors of all scenarios.These results show that the HAPS shows the best performance in terms of both estimated probability and localization accuracy.It is owing to the fact that the HAPS has a higher LoS ratio to guarantee a higher estimated probability, and also direct waves are more likely to arrive at the HAPS rather than multiple reflection paths that increase the angular estimation accuracy.

Conclusions
In this study, we compared the HAPS, UAV, and DEURAS systems for the location estimation of radio sources, and found that the HAPS estimation showed excellent accuracy at distances up to 20 km due to its high LoS ratio, while the UAV estimation showed poor performance due to interference from reflected waves, especially in urban environments.The position estimation using an antenna configuration that mimics the DEURAS system also showed lower position estimation accuracy compared to the HAPS due to the antenna height and reflected waves.Overall, the HAPS is excellent for monitoring radio sources over a wide area, and we believe that integrating the HAPS into the existing DEURAS setup has the potential to enhance the current system.
However, several problems are expected to arise when the HAPS is actually used for location estimation in practice.The first is the coordination among three HAPS units.The

Conclusions
In this study, we compared the HAPS, UAV, and DEURAS systems for the location estimation of radio sources, and found that the HAPS estimation showed excellent accuracy at distances up to 20 km due to its high LoS ratio, while the UAV estimation showed poor performance due to interference from reflected waves, especially in urban environments.The position estimation using an antenna configuration that mimics the DEURAS system also showed lower position estimation accuracy compared to the HAPS due to the antenna height and reflected waves.Overall, the HAPS is excellent for monitoring radio sources over a wide area, and we believe that integrating the HAPS into the existing DEURAS setup has the potential to enhance the current system.
However, several problems are expected to arise when the HAPS is actually used for location estimation in practice.The first is the coordination among three HAPS units.The DEURAS system uses a wide-area cabled LAN to aggregate information from each receiving antenna at the center station, but it is considered difficult to use this method for the HAPS since wired backhauling is impossible for aerial stations.However, such issues can be solved via wireless or optical backhauling among HAPS units, or by using feeder links to coordinate with other HAPSs via ground antennas.The second is the case where there are multiple sources of emitters in the same service area.As a first step to demonstrate the benefit of the proposed system, this paper assumed only a single source.But in realistic circumstances, there could be multiple sources emitted at the same frequency and time due to the nature of illegal emitters.In such a case, the HAPS units need to identify multiple sources simultaneously.This issue can also be solved via sharing the antenna resources equipped on the HAPS units.The third is when the illegal emitter is not fixed is but moving like a car or a truck.In such case, the introduction of a tracking filter upon the mechanism in this paper might be a solution.Although we have not yet been able to study these issues in depth, they will all be considered for our future works.
Furthermore, we used MATLAB to simulate the location estimation of a radio source.It is possible that diffracted and scattered waves, in addition to other arriving and reflected waves, may exist in a real environment, and these may affect the estimation results.Therefore, we would like to conduct actual experiments in the future to confirm the accuracy of the estimation.In addition, we would like to consider the location estimation of low-power radio sources, moving radio sources, and indoor radio sources, which were not covered in this study.

Figure 1 .
Figure 1.Number of wireless stations in Japan.

Figure 1 .
Figure 1.Number of wireless stations in Japan.

Figure 2 .
Figure 2. Statistics of illegal radio stations in Japan.

Figure 5 .
Figure 5.An example of an antenna pattern.

Figure 5 .
Figure 5.An example of an antenna pattern.

Figure 5 .
Figure 5.An example of an antenna pattern.

Figure 8 .
Figure 8. Simulation overview of LoS ratio calculation.

100mFigure 7 .Figure 7 .
Figure 7. 3D model of the evaluation environment.An overall view of the simulation is shown in Figure8.The simulation was conducted assuming that the HAPS was 2000 m above the ground and the UAV was 200 m above the ground.

Figure 8 .
Figure 8. Simulation overview of LoS ratio calculation.

Figure 10 .
Figure 10.An example of ray tracing simulation in progress.

Figure 10 .
Figure 10.An example of ray tracing simulation in progress.

Figure 10 .
Figure 10.An example of ray tracing simulation in progress.

Figure 11 .
Figure 11.Analysis of LoS ratio when sensors are away from the evaluation center.

Figure 11 .
Figure 11.Analysis of LoS ratio when sensors are away from the evaluation center.

Figure 12 .
Figure 12.Relationship between distance and LoS ratio.

Figure 13 .
Figure 13.Simulation overview of our location estimation process.

Figure 13 .
Figure 13.Simulation overview of our location estimation process.

Figure 14 .
Figure 14.Cylinder antenna used for simulation.

Figure 14 .
Figure 14.Cylinder antenna used for simulation.

Figure 15 .
Figure 15.Direction of arrival estimation using high altitude sensors.

Sensors 2024 , 24 Figure 17 .
Figure 17.CDF of errors using the HAPS at different distances against the emitting sources.

Figure 17 .
Figure 17.CDF of errors using the HAPS at different distances against the emitting sources.

Figure 18 .
Figure 18.CDF of errors using UAV at each distance.

Figure 18 .
Figure 18.CDF of errors using UAV at each distance.

Figure 20 .
Figure 20.CDF of errors for different sensor altitudes with sensor radius of 500 m.

Figure 20 .
Figure 20.CDF of errors for different sensor altitudes with sensor radius of 500 m.
E[•] denotes an ensemble average operation, λ is the wavelength of the carrier wave, d array is the antenna spacing, Re{•} donotes the real part, and Im{•} denotes imaginary

Table 2 .
Specifications of the cylinder antenna.

Table 3 .
Results with the HAPS 5 km away.

Table 4 .
Results with the HAPS 10 km away.

Table 5 .
Results with the HAPS 20 km away.

Table 6 .
Results with UAV 5 km away.

Table 7 .
Results with UAV 10 km away.

Table 6 .
Results with UAV 5 km away.

Table 7 .
Results with UAV 10 km away.

Table 8 .
Results with UAV 20 km away.

Table 8 .
Results with UAV 20 km away.

Table 9 .
Results with antenna placement assuming the DEURAS system.

Table 10 .
Results with the HAPS 500m away.

Table 9 .
Results with antenna placement assuming the DEURAS system.

Table 10 .
Results with the HAPS 500 m away.

Table 11 .
Results with UAV 500 m away.